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Abstract 


We  present  a  new  adaptive  genuinely  multidimensional  method  within  the  frame¬ 
work  of  the  discontinuous  Galerkin  method.  The  discontinuous  evolution  Galerkin 
(DEG)  method  couples  a  discontinuous  Galerkin  formulation  with  approximate  evo¬ 
lution  operators.  The  latter  are  constructed  using  the  bicharacteristics  of  multidi¬ 
mensional  hyperbolic  systems,  such  that  all  of  the  infinitely  many  directions  of  wave 
propagation  are  considered  explicitly.  In  order  to  take  into  account  multiscale  phe¬ 
nomena  that  typically  appear  in  atmospheric  flows  nonlinear  fluxes  are  split  into  a 
linear  part  governing  the  acoustic  and  gravitational  waves  and  to  the  rest  nonlinear 
part  that  models  advection.  Time  integration  is  realized  by  the  IMEX  type  ap¬ 
proximation  using  the  semi-implicit  second-order  backward  differentiation  formulas 
(BDF2)  scheme.  Moreover  in  order  to  approximate  efficiently  small  scale  phenom¬ 
ena  adaptive  mesh  refinement  using  the  space  filling  curves  via  AMATOS  function 
library  is  applied.  Three  standard  meteorological  test  cases  are  used  to  validate 
the  new  discontinuous  evolution  Galerkin  method  for  dry  atmospheric  convection. 
Comparisons  with  the  standard  one- dimensional  approximate  Riemann  solver  used 
for  the  flux  integration  demonstrate  better  stability,  accuracy  as  well  as  reliability 
of  the  new  multidimensional  DEG  method. 
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1  Introduction  and  Meteorological  Motivation 


A  characteristic  property  of  atmospheric  flows  is  their  multiscale  nature  with  wave  speeds 
differing  by  orders  of  magnitude.  If  the  Mach  and  Froude  numbers  are  small,  the  acoustic 
and  gravitational  waves  are  much  faster  than  advection,  but  only  the  latter  is  of  primary 
interest  for  numerical  weather  prediction.  Naive  explicit  time  integration  would  yield 
prohibitively  expensive  numerical  simulations,  which  makes  a  suitable  splitting  of  fast  and 
slow  waves  highly  desirable.  This  idea  is  not  new  and  has  been  used  extensively  in  previous 
meteorological  simulations.  Many  operational  nonhydrostatic  weather  models  use  split- 
explicit  methods,  where  different  time  steps  are  used  for  slow  and  fast  waves,  respectively, 
cf.  [14],  [23],  the  National  Center  for  Atmospheric  Research  [22],  Pennsylvania  State 
University /National  Center  for  Atmospheric  Research  [41]  or  German  Weather  Service 
[37].  Another  common  approach  is  based  on  semi- implicit  time  discretization;  here  the 
fast  waves  that  are  of  less  interest  are  approximated  implicitly,  whereas  slow  advection  is 
treated  explicitly.  Several  methods  following  this  idea  can  be  found,  e.g.,  in  [2,  8,  14,  20, 
21,  31,  30,  35,  40]  to  name  just  a  few. 

Another  characteristic  of  many  atmospheric  flows  is  their  multidimensional  character  with 
different  localized  structural  phenomena  such  as,  e.g.,  the  cloud-environment  interface.  A 
convenient  tool  to  approximate  these  local  structures  efficiently  is  mesh  adaptivity.  In¬ 
deed,  adaptive  mesh  refinement  has  been  applied  in  the  atmospheric  sciences  quite  suc¬ 
cessfully  over  the  past  decades,  see,  e.g.  [38,  3,  5].  Of  course,  the  strategy  where  and  how 
the  mesh  has  to  be  refined  is  a  difficult  scientific  problem  and  depends  on  the  particu¬ 
lar  application.  The  final  application  we  have  in  mind  is  the  simulation  of  an  evolving 
cumulus  cloud  and  its  interaction  with  the  environment.  This  is  an  important  meteo¬ 
rological  problem,  since  the  evaporative  cooling  at  the  cloud-environment  boundary  and 
its  impact  on  the  cloud  evolution  are  not  well  understood  [30],  [18].  Consequently,  effi¬ 
cient  adaptive  numerical  schemes  can  be  expected  to  improve  the  insight  into  underlying 
physical  processes  by  explicitly  resolving  the  interplay  between  the  larger  scales  of  the 
cloud  environment  and  the  smaller  scales  inside  the  cloud  and  at  its  boundary.  In  order 
to  approximate  localized  structures  efficiently,  we  will  work  with  adaptive  meshes  using 
the  space  filling  curves  via  the  AMATOS  function  library,  cf.  [4], 

In  this  paper  we  develop  a  new  semi-implicit  genuinely  multidimensional  method  within 
the  framework  of  discontinuous  Galerkin  scheme.  The  method  is  implemented  in  the 
discontinuous  Galerkin  solver  by  Giraldo  and  Warburton  [15],  see  also  recent  results 
[29,  30]  for  applications  to  the  Euler  equations.  However,  instead  of  a  standard  one¬ 
dimensional  approximate  Riemann  solver,  the  flux  integration  within  the  discontinuous 
Galerkin  method  is  now  realized  by  means  of  a  genuinely  multidimensional  evolution  op¬ 
erator.  The  latter  is  constructed  using  the  theory  of  bicharacteristics  in  order  to  take  all 
infinitely  many  directions  of  wave  propagation  into  account.  The  approximate  evolution 
operator  can  be  interpreted  as  a  multidimensional  numerical  flux  function.  In  the  finite 
volume  framework  the  finite  volume  evolution  Galerkin  (FVEG)  method  has  been  used 
successfully  for  various  physical  applications,  e.g.,  wave  propagation  in  heterogeneous 
media  [1],  the  Euler  equations  of  gas  dynamics  [6,  27]  or  the  shallow  water  equations 
[10,  19,  24],  The  FVEG  method  has  been  shown  to  be  more  accurate  than  standard  FV 
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methods  based  on  one-dimensional  Riemann  solver,  see  also,  [25]  for  further  references. 
One  interesting  question  arising  in  this  context  is  the  following:  does  the  multidimensional 
evolution  operator  bring  some  advantages  also  in  the  discontinuous  Galerkin  framework? 
In  order  to  illustrate  high  accuracy,  stability  and  robustness  of  the  new  DEG  method  we 
will  concentrate  on  two-dimensional  dry  atmospheric  flows  and  standard  meteorological 
test  cases. 


The  paper  is  organized  as  follows.  In  the  next  section  we  will  describe  the  mathematical 
model  governing  dry  atmospheric  flow,  in  Section  3  we  derive  the  discontinuous  Galerkin 
method  for  spatial  discretization  and  the  IMEX  type  approximation  for  time  discretiza¬ 
tion.  An  emphasis  is  put  on  a  non-standard  discretization  of  the  cell  interface  fluxes 
by  means  of  multidimensional  EG  operators,  cf.  Subsection  3.2.  Numerical  experiments 
presented  in  Section  4  illustrate  high  accuracy  and  stability  of  the  new  EG  method  and 
show  comparisons  with  the  DG  method  that  uses  the  standard  Rusanov  numerical  flux. 


2  Mathematical  Model 


We  start  with  the  description  of  the  mathematical  model.  Motion  of  compressible  flows 
is  governed  by  the  Euler  equations 


dtp  +  V  •  (pu)  =  0 

dt{pu)  +  V  •  (pu<g)  u  +  pld)  =  —  ppk  (2.1) 

dt{pO )  +  V  •  (pdu)  =  0  , 

where  p  denotes  the  density,  u  velocity,  p  pressure  and  9  the  potential  temperature. 
Further,  g  represents  the  gravitational  constant,  Id  is  the  identity  matrix  and  k  the  unit 
vector  in  the  vertical  direction.  Denoting  T  temperature,  the  potential  temperature  can 
be  obtained  from  the  equation  of  adiabatic  process  in  an  ideal  gas 

9  =  T  ^  /P,  R  =  cp-cv. 


We  use  potential  temperature  as  a  variable  since  it  is  better  suited  for  generalization  to 
moist  atmospheric  flow.  In  order  to  close  the  system  we  determine  pressure  from  the  state 
equation 

fRpoy 


P  =  P  0 


V  Po  J 


where  7  =  cp/cv  is  the  adiabatic  constant  and  p0  =  10 5 Pa  the  reference  pressure. 


Many  geophysical  flows  can  be  considered  as  a  perturbation  of  some  reference  equilibrium 
state.  For  example,  atmospheric  flows  are  typically  represented  as  a  perturbation  over 
the  background  hydrostatic  state  (p,  u(=  0),p,9),  cf.,  e.g.  [14],  [30], 


dp 

dy 


-pg- 
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Here  we  assume  that  9  =  9(y)  or  a  particular  case  9  =  300K.  Then  p  =  r  ,p  = 
p(p,  9)  =  p0  with  the  Exner  pressure  7 t(y)  :=  1  —  gy/cp9. 

In  order  to  avoid  numerical  instabilities  due  to  the  multiscale  behaviour  of  (2.1)  numerical 
simulations  are  typically  realized  for  perturbations  p'  =  p  —  p,  9'  =  9  —  9,  p'  =  p  —  p.  The 
latter  satisfy  the  following  equation 

dtp'  +  V  •  (pu)  =  0 

dt(pn)  +  V  •  (pu  ®  u  +  p'Id)  =  —  p'gk  (2.2) 

dt{p9)'  +  V  •  (p9u)  =  0. 

Our  aim  in  what  follows  is  to  approximate  (2.2)  with  the  discontinuous  Galerkin  method. 
However,  instead  of  the  classical  one-dimensional  numerical  flux  function  we  will  apply 
a  genuinely  multidimensional  evolution  operator.  To  this  end  let  us  rewrite  (2.2)  in  the 
form  of  hyperbolic  balance  law 


dq 

~dt 


+  V  •  F(q) 


S(q), 


(2.3) 


where 


and 


°  \ 
-p'<?k 

0  J 


(  PU  \  ( 

F(q)  =  pu®u  +  p'Id  ,  S(q)  = 

\  P9n  )  V 


is  the  nonlinear  flux  function  and  the  source  term,  respectively.  We  should  note  that  in 
our  numerical  experiments  we  will  also  use  a  stabilization  through  the  artificial  viscosity 
[39],  [30],  which  results  in  the  following  source  term 


S(q) 


0  \ 

-p'pk  + V  ■  (ppVu)  , 

V  •  (ppW)  J 


where  p  >  0  is  an  artificial  viscosity  parameter  that  will  be  specified  later. 


3  Discontinuous  Galerkin  method  and  the  multidi¬ 
mensional  EG  operator 

In  the  last  decades  the  discontinuous  Galerkin  (DG)  method  has  been  used  extensively  for 
approximation  of  partial  differential  equations,  see,  e.g.,  [7,  9,  11,  17,  32,  36]  and  the  ref¬ 
erences  therein.  The  method  is  popular  clue  to  its  flexibility:  it  is  based  on  the  Galerkin, 
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i.e.  variational,  formulation,  allows  handling  unstructured  triangulations  of  complex  ge¬ 
ometries,  mesh  adaptation,  and  flexible  choice  of  the  polynomial  basis.  Consequently, 
both  mesh  (h—)  and  polynomial  (p—)  adaptivity  can  be  applied  quite  naturally. 

In  this  section  we  follow  [35],  [29],  [30]  and  derive  the  strong  formulation  of  (2.3).  Let  us 
divide  the  computational  domain  12  into  a  finite  number  of  mesh  cells  Vte  with  a  boundary 
d Vte.  In  our  numerical  experiments  we  work  with  triangular  elements  12e  and  use  the  nodal 
basis  functions  {/</y ,  j  —  1, . . . ,  N},  N  is  a  number  of  degrees  of  freedom.  Now,  multiplying 
(2.3)  with  a  nodal  basis  ^(x),  integrating  over  12e  and  applying  twice  integration  by  parts 
we  obtain  the  strong  formulation 

(jit  +  V  '  F(q^  “  S(q^)  ^(x)dx  =  (3-4) 

/  (F(q,)-F*)^(x)d S,  i  =  1, . . . ,  N. 

Jdne 

Here  c\h  denotes  a  numerical  solution  q^(x)  :=  Q and  F*  is  a  suitable  numerical 
flux  function  that  approximates  cell  interface  fluxes. 

We  use  the  Lagrange  polynomials  having  the  basis  functions  ifjj  with  the  Fekete  points  for 
the  interpolation  and  the  Gauss  points  for  the  integrations.  In  the  numerical  experiments 
presented  below  approximations  with  the  second  order  Lagrange  polynomials  will  be  used. 
In  the  previous  papers  [13,  15,  29,  30]  ,  where  geophysical  flows  were  simulated  by  the 
discontinuous  Galerkin  method  the  one-dimensional  Rusanov  flux  has  been  chosen  as  the 
numerical  flux,  i.e. 

F*  :=  \  (F(qL)  +  F(qR))  -  A  (qR  -  qL)  n.  (3.5) 

Here  qL,  q^  denote  the  limiting  left  and  right  values  of  the  approximate  solution  at  the 
cell  interface,  n  is  the  outer  normal  to  the  cell  interface  and  A  the  maximum  wave  speed 
||u||  +  a,  a  denotes  the  speed  of  sound. 

The  novelty  of  our  work  relies  on  the  use  of  multidimensional  evolution  operator  in  order 
to  compute  F*.  More  precisely,  in  the  predictor  step  the  cell  interface  values  q*  are 
computed  by  the  multidimensional  evolution  Galerkin  operator  (EG),  q*  :=  EGqh.  The 
multidimensional  numerical  flux  is  then  obtained  as  follows 

F*:=F(q*).  (3.6) 

Derivation  of  the  multidimensional  evolution  operator  will  be  presented  in  Section  3.2. 


3.1  IMEX  time  integration 


In  order  to  approximate  efficiently  multiscale  phenomena  we  follow  the  ideas  of  Restelli 
and  Giraldo  [12],  [35]  and  split  the  system  (2.3)  into  two  parts  governing  fast  and  slow 
waves.  The  fast  waves  are  approximated  implicitly  and  the  slow  waves  explicitly.  More 
precisely, 


<9q 

~dt 


•A/"(q), 
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where  the  full  nonlinear  operator  A/"(q)  =  V  •  F(q)  —  S(q)  is  split  in  the  following  way 

A/"(q)  =  £(q)  +  K{  q) 


with  7 Z  :=  J\f  —  C  and 


C(q)  := 


(  div(ptt)  \ 
dp' /dx 
dp' /dy  +  gp' 

\  di  v(9pu)  ) 


Here  we  use  a  linearized  version  of  gradients  of  p'  and  set  ^  =  7xlz~’  where 

7  =  7-R  =  const.  Consequently,  the  operator  C  is  indeed  a  linear  operator  with  respect 
to  the  variables  q  :=  ( p',pu,pv ,  (p6)')T .  Denoting  S(q)  :=  (0,  0,  gp',  0)r  the  semi-discrete 
form  of  the  linear  subsystem  reads 


dq  T  dq  dq 

m+Jldi  +  J2dy 


S(q). 


Here  the  Jacobians  are 


(  o 

1 

0 

0  ^ 

0 

0 

0 

7 

0 

0 

0 

0 

^0 

9 

0 

o  ) 

(  0 

0 

1 

°  ^ 

0 

0 

0 

0 

0 

0 

0 

7 

0 

9 

o  ) 

(3.7) 


where  9  =  9(y). 

Now,  computing  the  eigenstrnctnre  of  the  matrix  pencil  of  the  system  (3.7)  shows  that 
(3.7)  is  a  hyperbolic  system  with  eigenvalues  Ai  =  — c,  A2,3  =  0,  A4  =  c,  where  c  :=  \f^9. 
It  should  be  pointed  out  that  for  the  non-dimensional  form  of  (2.2)  we  would  have  7  =  )g|, 
where  M  =  uref/cref  is  the  Mach  number  and  the  system  (3.7)  indeed  models  fast  acoustic 
waves,  see  also  [34],  Consequently,  it  has  to  be  approximated  in  an  implicit  way.  In  what 
follows  we  use  the  second  order  backward  difference  formula  of  implicit-explicit  (IMEX) 
type 


X  =  X  WK~m)  -  £(«;;■’”)]  +  a«;;+1).  p.s) 

m=—l  m= 0 

where  for  the  fixed  time  step  At  we  have  a_i  =  —  l,ao  =  4/3,  a\  =  —1/3,7  =  2/3,  (3o  = 

2,A  =  -l- 

The  scheme  can  be  also  rewritten  as  the  explicit  predictor  and  implicit  corrector  scheme, 
which  yields 


i  l 

9?  ~  X  a'”«r’n  +  X  (3-9) 

m=0  m=0 

and 
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(3.10) 


1 

[1  -  7A tC]  ql*'  =  q?  -  "I At  Y  /3m£(,rm). 

m= 0 

The  resulting  system  of  linear  algebraic  equations  is  solved  by  a  suitable  algebraic  solver. 
In  our  experiments  presented  below  we  have  used  the  GMRES  with  the  Jacobi  precon¬ 
ditioner,  more  sophisticated  choices  of  linear  solvers  may  further  decrease  computational 
costs. 


3.2  Multidimensional  EG  operators 


The  evolution  Galerkin  method  has  been  firstly  proposed  by  Lukacova,  Morton  and  War- 
necke  [26]  for  linear  acoustic  equation  and  later  generalized  by  Lukacova  and  coworkers 
in  the  framework  of  finite  volume  evolution  Galerkin  schemes  for  more  complex  hyper¬ 
bolic  conservation  laws,  such  as  the  Euler  equations  of  gas  dynamics  [27],  the  shallow 
water  equations  [24,  19,  10]  just  to  mention  a  few  of  them.  Since  the  flux  integrals  are 
approximated  using  the  multidimensional  evolution  operators  the  interaction  of  complex 
multidimensional  waves  is  approximated  more  accurately  in  comparison  to  schemes  us¬ 
ing  just  one-dimensional  approximate  Ricmann  solvers.  Extensive  numerical  experiments 
also  confirm  good  stability  as  well  as  high  accuracy  of  the  evolution  Galerkin  methods 
[26,  27,  24,  19], 

In  this  subsection  we  will  derive  approximate  evolution  operators  for  both  operators  A f 
as  well  as  C  that  are  based  on  the  theory  of  bicharacteristics  for  multidimensional  hyper¬ 
bolic  conservation  laws.  We  will  describe  the  derivation  of  the  evolution  operator  for  the 
operator  J\f  in  more  details,  the  derivation  of  the  evolution  operator  C  is  analogous. 

First,  let  us  rewrite  (2.2)  in  a  quasilincar  form  using  the  primitive  variables  w  =  (p',  u,  v,p') 


dtw  +  A^w)  dxw  +  A2(w)  cLw  =  s(w) 


(3.11) 


with 


A1  = 


P  0  0  \ 

o 

CL 

o 

(  °yp  V  \ 

u  0  -p 

o 

o 

o 

0 

0  u  0 

,  A-2  — 

o 

o 

ce 

It- 

,  S  =  - 

—g 

p 

7P  0  u  ) 

\  0  0  7  p  v  ) 

\  9yPV  ) 

(3.12) 


Using  the  above  thermodynamic  relationship  for  p,p  we  obtain 


dyp  = 


P  0  C-v  Q 

(Rd)2cp 


1  - 


gy 

Cp  6 


CV  _ 1 

R  1 


8yP  = 


gpo 

Re 


i  - 


gy 

Cp  9 


cV 

R 


We  first  linearize  (3.11)  by  freezing  the  Jacobian  matrices  Al7A2  at  a  suitable  interme¬ 
diate  state  p',  u,  v,  p' .  In  the  numerical  experiments  presented  below  this  intermediate 
state  is  obtained  as  the  local  average  of  (left /right)  limiting  values  at  cell  interfaces.  Since 
our  problem  is  hyperbolic  we  have  real  eigenvalues  and  a  full  set  of  linearly  independent 
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P  =  (x,y,tn  +  T) 


Figure  1:  Bicharacteristic  cone  used  for  the  evolution  operator 

eigenvectors  corresponding  to  the  matrix  pencil  P  :=  Ainx  +  A2ny,  where  (nx,ny)  is  an 
arbitrary  unit  vector  ||(nx,nl,)||  =  1.  Indeed,  the  eigenvalues  of  P  are 

Ai  =unx  +  v7iy  —  a , 

A2  =  A3  =  unx  +  vny, 

A4  =  unx  +  v  ny  +  a, 


where  a  :=  J7?  =  \/ 7 R6 


pR9  \ 

Po  J 


is  the  speed  of  sound.  Now  multiplying  (3.11)  by  a 


matrix  R  1,  R  consists  of  the  right  eigenvectors  of  P,  we  can  rewrite  (3.11)  using  the 
so-called  characteristic  variables  v  =  R 


w 


<9tv  +  Bxdxy  +  B2dy\  =  r  , 

where  Bi  :=  R  1  A[R  and  r  :=  R  1  s(w).  Equivalently,  we  have 

dty  +  diag(S1)93,v  +  dmg(B2)dyv  =  S  +  r  (3.13) 


with 

S(x,  9)  :=  -  [{B1  -  diag(Si))c)xv  +  (B2  -  diag(S2))^v] . 

Integrating  each  equation  of  (3.13)  along  the  corresponding  bicharacteristic  x3  from  tn  to 
tn  +  r: 

d*j 
dt 


[-® 1  ,n  5  -®2 ,jj]  1  j  1)  •  •  •  >  4, 


we  obtain  after  some  lengthy  manipulations,  analogous  to  those  in  [26,  27],  the  following 
exact  integral  representation 


p\  pi  = 


p 

2nd 


p2tt 


cos(cn)  w(Qi(a;))  —  sin  (a;)  n(Qi(a;))  +  — p'fQiku)) 

p  a 


du 


+  ✓(  <+)-PW- 

n  * 


f 3(t,u )  dtdu 

—  sinfadp— (xi(f, a;))  H — ^  ^  ’ — —dypdtdu 
p  pa 


rtn+T 

+L  v 


(x2(f))  dyp+  ^ J  dt, 


2tt  r  u 


“<P|  -  2n 


p'(  QiM) 

pa 


cos(cn)  +  w(Q1(a;))  cos2(cn)  +  u(Q1(a;))  sin  (a;)  cos(cn) 


den 


1  1 

“tl  (Q2)  +  7^ —  /  /  cos(cn)  /3(f,  u)  dt  du 

2  2?r  do  dtn 

1  /“2,T  /‘fn+r  u(xi(f,a;))  _  p' 

—  /  /  cos(cn) - — - OyP  —  sm(w)  cos  (a;)  <7  —  (xi(£,  enj)  at  acn 

2?r  do  dtn  pa  P 

2  r^nJrT 


—  J  dxp'(x2(t))dt , 


«(P)  =  FZ 


1 

2~7T 


2^  r 


p'iQi) 

pa 


sin(cn)  +  n(Q1)  cos(tn)  sin  (a;)  +  u(Qi)  sin2(cn) 


du 


Y  Y  /“27r  ptn+T 

-u(Q2)  +  —  /  /  sin  (u)  p(t,u)dtdu 


JO  J  t-, 

^  p2n  ftn+T 


arJn  ./,  sm(w'  PS 


v(xi(t,  en))  .  200..P' 


'0  Jt 
rtn+T 


dyp  —  siir(tn)p— (xx(f,  u))  dt  du 


—  jt  dyp'(x2{t))  dt  -  -g 


1  /^+T  p' 


—  (x2(t))dt, 
P 


r»27r 


p'(P)  =  —  /  tp'(QiH)  —  paw(Q1(a;))  cos(tn)  —  pau(Q1(a;))  sin(cn)]  du 
2n  J o 

+T 

/3(t,  u)  dt  du 

—  sm(u)pag—(xi(t,u))  +  v(x.i(t,u))dypdtdu. 


ffere  / 3(t,u )  =  a[dxu  sin2(cn)  —  (dyu  +  dxv)  sin(cn)  cos (u)  +  dyv  cos2(cn)]  and  P  =  (x,y,tn  + 
r),  Qx  =  Qi(cn)  =  (x  —  (u  —  a  cos(tn))r,  y  —  (v  —  dsin(cn))r,  tn),  Q2  =  (x  —  ut,  y  —  hr,  tn) 
are  respectively  the  pick  and  footpoints  of  the  bicharacteristics  that  generate  the  mantle 
of  the  bicharacteristic  cone,  cf.  Figure  1. 
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To  obtain  an  explicit  approximate  evolution  operator  the  above  exact  integral  represen¬ 
tation  needs  to  be  approximated.  First,  all  time  integrals  are  approximated  using  the 
rectangle  rule.  Integrals  along  the  base  perimeter,  that  contain  (3(tn,uj)  terms,  are  re¬ 
placed  by  means  of  integration  by  parts.  More  precisely,  applying  Lemma  2.1  [26]  we 
obtain 

/*27T  /*27T 

r  /  /3(tn,uj)  =  /  ^(Q-l)  cos(u;)  +  u(Qx)  sin  (a;)  du 

Jo  Jo 


r*27r 


-2tt 


r  /  /3(tn,  u)  cos(a;)  =  /  u(Q1)(2  cos2(cu)  —  1)  +  2u(Q1)  sin  (a;)  cos(a;)  du 


-2tt 


-2tt 


r  /  /3(tn:u)  sin  (a;)  =  /  2u(Q1)  sin  (a;)  cos(a;)  +  u(Q1)(2  sin2(a;)  —  1)  du. 


Further  we  approximate  ^  with  4  and  substitute  the  condition  for  hydrostatic  balance 
dyp  =  —pg  in  the  equations  for  p'(P),  w(P),  v(P),p'(P).  The  state  equation  also  yields 


dyp  =  -7r" 


c-v  -  pg 

c„  R6 ^  a 2 


This  implies  that  the  last  term  in  the  equations  for  p'(P)  vanishes  approximately 


v(x2(t))  f-dyp-h  )  dtra  0. 
\  CL 


Moreover,  pressure  derivatives  dxp'(x. 2)  and  dyp'{x2)  in  the  equations  for  u(P),u(P)  can 
be  approximated  using  the  midpoint  rule  and  the  Gauss  theorem 

r  1  f2n 

-77zdxp\ Q2)  =  —  /  p'(Qi)  cos(o;)  du  +  G>(r2) 

2 p  2npa  J0 

t  i  r2n 

~7pdyP\ Q2)  =  -7T~. ~  /  P'(Qi)  sinM  du}  + 

zp  Z7T pa  Jq 

These  approximations  yield  finally  the  desired  approximate  evolution  operator  in  order 
to  predict  the  cell  interface  values  q*  =  (//(P),w(P),v(P),p'(P))  =  EG^  qh  in  (3.4). 
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Approximate  evolution  operator  EG ‘v  for  the  nonlinear  operator  M  (3.11) 


7(P)  = 


r»27r 


2ti  a 


—2pcos(u)  ^(Q^a;))  —  2psin(o;)  ^(Q^a;))  + 


P'(QiM) 


du 


P  (Q2) 


P'( Q2)  ,  r 


r»27 r 


+ 


2tt  a 


sin(o;)p/(Q1)  g  +  pg 


du , 


u(P)  =  — 


1 

2tt 


2tt  r  r>  / 


2p/(Qi; 

pa 


cos(o;)  +  m(Q1)(3  cos2(o;)  —  1)  +  3v(Q1)  sin  (a;)  cos(a;) 


du 


1  r  f27T 

-2^-2 rd  1 


cos  u 


sin(a;)p/(Q1)^  +  l^-^-pg 


du , 


»<p)  =  ^ 


i  r 27r  r  0/7  fo  i 

sin(a;)  +  3m(Q1)  sin(a;)  cos(a;)  +  (3sin2(a;)  —  l)v(Q1) 


pa 


du 


T 


r»27r 


sin(w)7(Qi)  g  + 


T 


duj  ~  2pP'^9  ’ 


r»27r 


p'(P)  =  —  /  [p,(Q1(a;))  —  2pera(Q1(a;))  cos(a;)  —  2pav(Q1(u))  sin(a;)]  du 

Jo 


t  a 
~2n 


"27T  T  v(Ci  ) 

sin(a;)p/(Q1)  g  +  Lpg 


du. 


(3.14) 


We  should  point  out  that  all  integrals  along  the  base  perimeter  (sonic  circle),  i.e.  integrals 
with  respect  to  u,  are  evaluated  exactly  for  given  discrete  data.  We  make  a  transformation 
of  the  actual  triangle  to  the  reference  triangle,  where  the  corresponding  integrals  along 
the  arcs  of  sonic  circle  were  precomputed  with  the  help  of  the  computer  algebra  package 
Mathematica. 


Our  next  goal  is  to  derive  the  approximate  evolution  operator  for  the  system  (3.7).  Since 
the  acoustic  subsystem  (3.7)  is  linear,  the  derivation  of  the  exact  evolution  operator  is 
even  easier  and  it  has  been  already  obtained,  e.g.,  in  [27].  Indeed,  applying  the  above 
procedure  we  obtain  in  an  analogous  way  the  following  exact  evolution  operator  for  the 
acoustic  system  (3.7),  cf.  [27],  [28]. 
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p'(p)  =  p'(Q 


7(p^)'(Qs 


2/  ~2 
Cz 


r'2n 


±1  -  cos(w)  (pw)(Q1(w)  -  sin(a;)  (puXQiM)  +  ?(p£)/(QiM)dw 

JO  c 

2  ^»27T  l*tn+T  2 

[cos(a;)(pM)(xi(t, a;))  +  sin(a;)(/w)(xi(t, a;))]  dtdcu 


27TC 


2itc 


'0  Jtn  Lri 

/*27T  ptn-\-T 


tn  +  T-t 


sm(u)gp'(x.i(t,u))  dtdu , 


'o  Jt„ 


rtn+T 


(pu)(p)  =  ^(pm)(Q2)  -  |  /  ^(pP)'(x2(t))dt 

J  tn 

1  r2?r 

+  —  /  cos2(u;)  (pw)(Q1(a;))  +  sin  (a;)  cos(u;)  (pv)(Q1(a;))  —  -(p0)/(Q1(a;))  cos(u;)<icu 
2?r  Jo  c 

2  /*27T  rtn-\~T  ^ 

- /  /  - [cos(2a;)  (pw)(xi(f,  a;))  +  sin(2a;)  (pu)(xi(t,  a;))]  dt  du 

2vrc  y0  Jtn  tn  +  T  -t 

1  ^277  rtn-\~T 

H - /  /  sin  (a;)  cos(a;)pp/(xi(t,  a;))  dt  du  , 

2?rc  Jo  Jtn 


fin+T 


(H(P)  =  ^M(Q2)  -  |  /  c>J/(pP),(x2(t))dt 


-2tt 


—  /  cos(a;)  sin  (a;)  (pw)(Q1(a;))  +  sin2(a;)  (pv)(Q1(a;))  —  ^(pP),(Q1(a;))  sin(a;)ciu 
2tt  Jo  c 


rtn+T 


+ 


2vrc  ,/o  ytn  tn  +  T-t 

1  1‘2,'K  rtn+T 


sin(2a;)  (pw)(x!(t,  ca))  —  cos(2ca)  (pi;)(x!(£,  a;))]  dtdca 


27TC 


siir(a;)pp'(x1(t,  a;))  dt  du  , 


0  J  tn 


-2?r 


(pP)'(P)  =  /  (pP)/(Q1(P))  -  5  cos(ca)  (pu)(Q1(a;))  -  ?  sin(ca)  (pu)(Q1(w))dc1; 

2vr  J0  7  7 


p27T  f*tn~\~T 


2vr  Jo  Jtn  lf(jn  +  T-t ) 

2  p2?r  rtn+T 


cos(ca)  (p-u)(x1(t,a;))  +  sin  (a;)  (pi;)(x1(i,  a;))]  dt  du 


27T7 


sin(a;)pp,(x1(t,  u))  dt  du. 


’0  j  tn 


Here  c=  y  7#  and  0  is  obtained  as  a  local  average  of  9.  Note  that  we  have  used  analogous 
notation  for  bicharacteristics  xi(t,  u),  x2(t)  and  Ql5  Q2  as  for  the  nonlinear  operator 
above.  Of  course,  in  the  case  of  the  linear  operator  £  the  corresponding  eigenvalues 
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determine  the  geometry  of  the  bicharacteristic  cone.  Consequently,  we  now  have  Qx  = 
(x  —  ccos(u)r,  y  —  csin(o;)r,  tn),  Q2  =  ( x ,  y,  tn)  and  P  =  ( x ,  y ,  tn  +  r). 

In  order  to  approximate  the  above  exact  integral  representation  we  apply  the  rectangle 
rule  for  time  integrations  and  the  second  and  third  equation  of  (3.7)  to  eliminate  space 
derivatives  of  (p0)'(Q2(f))  in  the  equations  for  pu  and  pv.  This  yields  the  approximate 
evolution  operator  EGC  that  has  been  denoted  as  the  EG3  in  [28] .  Consequently,  we  have 
for  the  acoustic  system  (3.7)  q*  =  (p'(P),  pw(P),  pv(P),  (p0)'( P))  =  EGC  qh. 


Approximate  evolution  operator  EGC  for  the  linear  operator  £  (3.7) 

7(/^),(Q2) 


p'(P)  =  p'(Q2)  - 


1 


-2tt 


27t  a 
T 

2nd 


-2cos(ta)  (pu)(Q1(u)  —  2sin(o;)  (pn)(Q1(o;))  +  ^(pd),(Q1(ta)) 


du 


r»27r 


sin(o;)pp,(xi(t,  u))  du 


(pu)(P)  = 


1 
7 r 


T 

2nd 

p  2n 

'o 


r»27r 


sin  (a;)  cos(a;)pp'(x1(t,  u))  du 

(3cos2(a;)  —  1)  (pw)(Q1(a;))  +  3sin(a;)  cos  (a;)  (pn)(Q1(a;)) 

-  ?(p0)'(Q1M)cos(u;) 
a 


du 


0*0  (P)  = 


T 

2nd 

p2n 


f*2n 


sin2 (ca)pp/(x1(t,  u))  du 


n 


3  cos(ta)  sin  (a;)  (pu)(Q1(u))  +  (3sin2(u;)  —  1)  (pn)(Q1(a;)) 

-  ?(p0)/(QiM)sin(w) 
a 


du 


(P0)'( P)  = 


T 


d2ir 


l 

2 n 


2n 7 

f*  2n 


sin(a;)pp'(xi(t,  u))  du 


>o 


{p0)\ Qi(0))  -  2?  cos(tn)  (pn)(Q1(w))  -  2?  sin(ca)  (pu)(Qi(w)) 
7  7 


du. 

(3.15) 


The  resulting  discontinuous  evolution  Galerkin  scheme  is  based  on  the  spatial  discretiza¬ 
tion  (3.4),  IMEX  type  time  discretization  (3.9),  (3.10)  and  the  multidimensional  approx¬ 
imate  evolution  operators  (3.14)  and  (3.15)  in  order  to  approximate  cell  interface  fluxes 
for  the  nonlinear  and  linear  operators  J\f  and  £,  respectively.  The  parameter  r  describes 
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the  time  step  of  the  evolution.  In  our  numerical  experiments  we  have  taken  r  <C  At, 
where  At  is  the  time  step  of  the  IMEX  time  integration,  cf.  Section  3.1.  Consequently, 
the  operator  models  the  local  evolution.  In  the  numerical  experiments  presented  below 
piecewise  quadratic  approximate  functions  in  space  are  used,  that  yield  six  degrees  of 
freedom  (DoF)  per  a  triangular  finite  element:  three  DoFs  are  located  at  the  vertices  of 
the  triangular  mesh  cells  and  three  DoFs  at  the  centers  of  the  cell  interfaces.  When  we 
benchmark  our  model  comparing  it  to  the  DG  method  with  the  Rusanov  flux  function,  we 
do  not  use  any  kind  of  limiter  or  filter  except  for  the  artificial  viscosity  as  described  in  Sec¬ 
tion  2.  In  our  recent  work  [6]  we  have  also  generalized  the  discontinuous  Galerkin  solver 
of  Giraldo  and  Warburton  [15]  by  including  the  GPU  implementation  of  the  genuinely 
multidimensional  EG  operator. 


4  Numerical  Experiments 

To  verify  the  new  discontinuous  evolution  Galerkin  method,  we  carry  out  three  numerical 
experiments  for  the  test  cases:  (i)  free  convection  of  a  smooth  warm  air  bubble,  (ii) 
collision  of  a  large  warm  bubble  with  a  small  cold  bubble  placed  on  top  of  the  warm  one, 
and  (iii)  a  density  current  caused  by  a  falling  cold  air  bubble.  In  these  experiments,  clue 
to  the  buoyancy  forces  caused  by  differences  in  the  density  of  the  air  bubbles  and  the 
isothermal  environment,  the  initially  resting  air  bubbles  develop  a  vertical  motion  with 
low  Mach  number  M  <  1. 


4.1  Test  1:  Free  convection  of  a  smooth  warm  air  bubble  due 
to  Giraldo  and  Restelli  [13] 

In  the  Giraldo-Restelli  experiment  shown  in  Figure  2,  the  warm  bubble  rises  and  deforms 
symmetrically  due  to  the  shear  friction  with  the  surrounding  air  at  the  warm/cold  air 
interface,  adapting  a  mushroom-like  shape  gradually.  The  warm  air  bubble  is  placed  at 
xc  =  500m,  yc  =  350m  with  the  initial  potential  temperature  perturbation: 


O' 


0  for  r  >  rc 

{0'c/ 2)  [1  +  cos  (7rr/rc)]  for  r  <  rc 


(4.16) 


where  6'c  =  0.5  K  is  the  maximal  initial  amplitude  of  the  perturbation,  rc  =  250m  is  the 
bubble  radius,  and  r  the  distance  to  the  center  of  the  bubble  ( xc,yc ). 

In  order  to  simulate  efficiently  localized  structures  arising  in  geophysical  flows,  such  as 
cloud  boundaries,  adaptive  mesh  refinement  is  a  very  suitable  tool.  We  work  with  the 
function  library  AMATOS  of  Behrens  et  al.  [4],  where  h— adaptive  mesh  refinement  is 
based  on  the  space  filling  curve  approach.  Our  numerical  experiments  were  performed 
on  the  domain  of  1  km  x  1  km  with  no-flux  boundary  conditions,  using  the  h— adaptive 
mesh  refinement  method,  where  the  spatial  resolution  is  adapted  by  refining  or  coarsen¬ 
ing  the  mesh  cells.  The  maximal  resolution  degree  of  the  mesh  is  n  =  11,  which  yields 
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1000m/v/2,!+1  ~  15.6m  per  the  finite  element  shortest  edge  in  the  simulation  domain. 
Since  we  use  the  second  order  polynomials,  the  finest  resolution  of  the  grid  is  7.8m  corre¬ 
sponding  to  a  half  of  the  length  of  a  shortest  edge. 

Analogously  as  in  [30],  in  the  numerical  experiment  presented  in  Figure  2  we  use  a  slightly 
modified,  simple  refinement  criterion 

max  [sgn(0')0'(x,  t)]  >  a\6'c\  (4.17) 

for  the  deviation  of  the  potential  temperature  from  the  background  state  O'  =  9  —  0; 
a  <C  1  is  a  test  dependent  parameter  (for  the  numerical  experiments  in  this  work  we  use 
a  =  0.1),  and  9'c  is  the  maximal  initial  amplitude  for  the  perturbation  of  the  potential 
temperature. 

If  condition  (4.17)  holds  on  some  element  h2e,  the  element  will  be  recursively  refined  up 
to  a  specified  finest  mesh  resolution.  In  the  rest  of  the  computational  domain  the  mesh 
is  adaptively  coarsened,  see  also  [30]  for  further  details. 

In  Figure  2  we  show  the  simulation  results  obtained  using  the  newly  derived  multidi¬ 
mensional  DEG  method  (left-hand  side  of  the  shown  snapshots)  and  the  DG  method 
with  the  standard  Rusanov  numerical  flux  (3.5)  (right-hand  side).  The  latter  model  will 
be  simply  referred  to  as  the  Rusanov  method  later  on  in  this  paper.  Both  simulations 
were  performed  using  the  IMEX  type  approximation  for  time  discretization  as  described 
in  Section  3.1.  The  artificial  viscosity  in  both  tests  is  kept  constant  n  —  0.1  over  the 
computational  domain  and  time.  The  integration  time  step  is  At  =  0.1s.  The  corre¬ 
sponding  Courant-Friedrichs-Lewy  conditions  calculated  on  the  smallest  mesh  element  is 
CFL  =  max  ~  4.43  for  the  total  velocity,  CFLU  =  max  <  3.3  •  10~2  for 

the  advective  part.  The  time  step  for  the  evolution  of  EG  operator  is  r  =  1.25  •  10-3s, 
which  corresponds  to  CFLEG  a;  5.5  •  10~2. 

Both  models  yield  very  similar  solutions  of  the  Euler  equations  model  for  the  time  t  < 
600s.  The  differences  between  both  solutions  become  however  noticeable  at  t  =  600s. 
The  solution  obtained  with  the  Rusanov  flux  exhibits  short  length  scale  oscillations  (cf. 
isolines  at  x  ~  350m,  y  800  in  Figure  2,  right-hand  side,  t  =  600s)  whereas  the  isolines 
are  quite  smooth  for  the  DEG  method  (on  the  left-hand  side  in  the  same  figure).  A  very 
pronounced  oscillation  of  the  perturbation  front  can  be  observed  on  the  top  of  the  air 
bubble  for  the  Rusanov  flux.  The  snapshots  for  the  time  t  =  900s  in  the  same  figure  show 
that  the  solutions  become  very  different  for  later  time. 

In  order  to  analyse  in  more  detail  the  evolution  of  the  solutions  at  later  times  we  plot 
a  few  selected  iso-levels  of  the  excess  potential  temperature  for  intermediate  times  in 
Figure  3.  The  isolines  shown  for  6'  =  0.05,  0.25,  and  0.4  K  represent  the  structure  of  the 
solution  close  to  the  background,  in  the  middle,  and  near  the  maximum  of  the  background 
perturbation,  respectively.  The  solution  based  on  the  EG  operator  can  be  characterized 
as  being  stable  against  the  short  length  scale  perturbations  (all  isolines  are  smooth).  A 
large  scale  oscillation  appears  at  the  bottom  of  the  air  bubble  interface  at  time  t  =  800s 
(x  ~  850m,  y  «  600m)  that  seems  to  be  due  to  the  Kelvin-Helmholtz  instability  caused 
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Figure  2:  Excess  potential  temperature  6'  for  the  rising  thermal  bubble  in  the  numerical 
experiment  proposed  by  Giraldo  and  Restelli  [13],  obtained  by  the  new  semi- implicit 
method  on  the  h— adaptive  grid  with  the  coarse/fine  grid  resolution  levels  n  —  1  —  11, 
respectively,  and  the  constant  artificial  viscosity  /i  =  0.1m2/s:  On  the  left-hand  side  is 
with  novel  multidimensional  EG  operator  in  the  numerical  flux  function,  on  the  right-hand 
side  with  the  Rusanov  flux  function.  The  real-world  domain  is  1  km  x  1  km  (only  a  half  of 
the  squared  computational  domain  is  shown  in  the  x-direction);  the  shortest  edge  of  the 
adaptive  mesh  elements  is  ~  15.6m,  the  spatial  resolution  ss  7.8m.  The  simulation  times 
are  as  indicated;  CFL  «  4.43,  advective  CFLU  G  {0,0.0238,0.0331,0.0286}  for  DEG 
method  and  CFLU  G  {0,  0.0238,  0.0332,  0.0313}  for  the  Rusanov  flux  model.  CFLEq  ~ 
0.055.  Contour  levels  correspond  to  9'  =  0.05,  0.15,  0.25,  0.35,  and  0.45  K. 
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Figure  3:  Giraldo-Restclli  test  for  the  DEG  method  a)  and  the  DG  method  with  the 
Rusanov  flux  b),  with  constant  viscosity  //  =  0.1m2/s  at  times  as  shown;  CFL  4.43, 
advective  CFLU  G  {0.0331,  0.0312,  0.0278,  0.0286}  for  DEG  model  (from  top  to  bottom) 
and  CFLU  G  {0.0332,0.0314,0.0279,0.0313}  for  the  Rusanov  flux.  Contour  levels  are  for 
6'  =  0.05,  0.25,  and  0.4  K.  The  structure  of  the  adaptive  meshes  is  shown  once  for  clarity. 
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Figure  4:  Giraldo-Restelli  test  for  the  DG  method  with  the  Rusanov  flux,  with 
constant  viscosity  fi  =  0.25m2/s  at  times  as  shown,  CFL  m  4.43,  CFLU  e 
{0.0306,0.0299,0.0269,0.0248}  (from  top  to  bottom).  Contour  levels  are  for  &  =  0.05, 
0.25,  0.4  K. 
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Figure  5:  Comparison  of  the  contour  level  9'  =  0.05  K  of  the  excess  potential  temperature 
in  the  Giraldo  and  Restelli  test  for  time  t  =  900s  using  the  DEG  method  with  constant 
viscosity  /i  =  {0.1,  0.18,  0.25}m2/s  (red,  green,  and  yellow  curves,  respectively),  and  the 
Rusanov  flux  model  with  constant  viscosity  p  =  0.25m2/s  (blue  curves).  The  coarse/fine 
mesh  resolution  is  n  —  1  —  11,  corresponding  to  the  spatial  resolution  of  fa  7.8m  along  the 
shortest  element  edge  (15.6m),  CFL  fa  4.43,  CFLU  fa  0.0286  (DEG)  and  CFLU  fa  0.0248 
(Rusanov). 
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by  the  shear  flow  interactions  along  the  interface  between  the  rising  bubble  with  the 
surrounding  air  at  rest.  These  long  wave-length  oscillations  of  the  interface  lead  to  the 
turbulent  structure  at  later  times  (e.g.,  t  =  900s  shown  in  the  same  figure),  typical 
for  the  Kelvin-Helmholtz  phenomenon.  This  has  been  studied  intensively  by  numerical 
experiments,  see,  e.g.,  [16,  30,  33];  note  that  the  exact  solution  to  this  problem  is  unknown. 

We  compare  results  obtained  by  the  DEG  method  with  those  obtained  by  the  DG  method 
with  the  Rusanov  flux  model  shown  on  the  right-hand  side  in  Figure  3  for  which  the  isolines 
exhibit  high  frequency  oscillations,  as  mentioned  above.  At  late  time  (t  >  700s)  island¬ 
like  tiny-scale  isolated  spots  (mini-bubbles)  develop  due  to  the  reconnections  between  the 
close  located,  oscillating  isolines.  Such  structure  can  be  seen  for  all  iso-levels  (not  only 
for  that  one  located  close  to  the  background  level  9'  —  0,  where  the  oscillations  of  the 
background  could  yield  such  artifacts). 

A  simple  increasing  of  the  artificial  viscosity  in  order  to  damp  the  oscillations  cannot  really 
be  considered  as  a  reliable  method  to  eliminate  this  issue.  Although  it  is  known  that  by 
increasing  the  viscosity  one  can  damp  the  oscillations,  we  show  that  the  final  solution  can 
change  on  large  scales  significantly,  too.  In  Figure  4  we  show  the  simulation  results  using 
the  Rusanov  flux  model  with  the  constant  artificial  viscosity  /x  =  0.25m2/s.  This  value  of 
the  viscosity  has  been  chosen  because  it  is  just  enough  to  make  the  model  stable  against 
the  short  length  scale  fluctuations  and  hence  to  remove  the  high  frequency  oscillations. 
However,  this  value  of  fi  is  already  so  high  that  the  obtained  solution  looks  very  different, 
if  compared  to  the  previously  discussed  results  for  /x  =  0.1m2/s:  i)  the  Kelvin-Helmholtz 
instability  along  the  outer  interface  has  been  completely  damped,  and  ii)  the  solution 
became  very  dissipative  which  can  be  deduced  from  the  shape  of  the  isolines  for  6'  =  0.4 
(the  red  curve  in  Figure  4).  This  isoline  shrunk  in  length  for  times  t  =  600s  and  t  =  700s 
and  it  has  completely  disappeared  for  later  times  that  means  a  significantly  decreased 
amplitude  of  the  perturbation,  if  compared  to  ji  =  0.1m2/s  in  Figure  3. 

In  Figure  5  we  show  a  comparison  of  the  solutions  obtained  by  the  DEG  model  for  several 
values  of  the  viscosity  parameter  yfx  =  0.1,  0.18,  0.25m2/s  and  by  the  Rusanov  flux  model 
with  the  viscosity  parameter  /x  =  0.25m2/s  chosen  for  the  stability  reason.  The  results  of 
both  models  agree  very  well  in  some  regions  of  the  computational  domain,  e.g.,  along  the 
interfaces  located  in  the  interior  of  the  air  bubble  where  the  shear  forces  are  negligible. 
The  outer  interface  calculated  by  the  DEG  method  shows  significant  difference  in  the 
details  of  the  simulated  structure  for  the  low  value  of  the  artificial  viscosity  parameter 
(yxx  =  0.1nr/s)  and  it  approaches  the  solution  of  the  Rusanov  model  as  the  viscosity 
parameter  increases.  We  would  like  to  point  out  that  the  DEG  solution  remains  stable 
for  different  choices  of  the  constant  viscosity  parameter,  whereas  for  the  Rusanov  model 
either  higher  constant  viscosity  or  adaptive  viscosity  (4.18)  has  to  be  applied  to  stabilize 
the  solution. 


20 


4.2  Test  2:  Collision  of  a  large  warm  bubble  with  a  small  cold 
bubble  due  to  Robert  [33] 

Robert  [33]  has  proposed  the  experiment  shown  in  Figure  6.  The  shape  of  the  rising  warm 
bubble  is  affected  by  the  small  cold  bubble,  which  slides  downwards  along  the  right-hand 
side  of  the  interface,  destroying  the  symmetry  of  the  warm  bubble.  Two  bubbles  are 
initially  placed  at  (xc,yc)  =  (500m,  300m)  and  (. xc,yc )  =  (560m,  640m),  for  the  warm 
and  the  cold  bubbles,  respectively  (cf.  Figure  6).  The  maximal  amplitudes  of  the  initial 
potential  temperature  perturbation  are  9'c  =  0.5  K  and  9'c  =  — 0.15K,  respectively.  The 
profiles  of  the  initial  perturbation  for  the  excess  potential  temperature  are  given  by  a 
Gaussian  distribution 


q,  _  f  9'c  for  r  <  rc 

\  9'cex p  [—  (r  —  rc)2  / 502]  for  r  >  rc 

with  a  flat  core  of  radius  rc  =  150m  for  the  warm  bubble  and  rc  =  0  for  the  cold  bubble. 

In  the  previous  tests  we  have  seen  that  one  can  obtain  a  stable  and  an  oscillating  solution 
on  an  adaptive  mesh,  depending  on  the  value  of  the  viscosity.  A  high  viscosity  stabilizes 
the  solution,  but  changes  it  too  much.  Here  we  perform  the  simulations  with  both:  on  a 
regular  static  and  on  an  adaptive  mesh,  with  a  constant  and  an  adaptive  viscosity  in  order 
to  study  the  impact  of  these  approaches  on  the  stability  of  the  methods.  The  integration 
time  step  is  set  to  At  =  0.1s. 

We  start  with  the  static  regular  grid  of  the  resolution  level  n  —  8  and  an  adaptive  artificial 
viscosity  in  the  source  term,  as  mentioned  in  Section  3.  We  follow  the  approach  of  [30], 
[39]  to  stabilize  the  simulations  through  the  adaptive  artificial  viscosity  calculated  from 

pe  =  max  ^o,Airef^^72(12_n)/2^  ,  (4.18) 

where  po  is  the  constant  viscosity  parameter  given  by  the  test  case,  pref  =  (l/17.7)m2/s 
and  a  =  0.4  are  the  test  independent  empirical  parameters,  A9'e  is  the  maximal  deviation 
of  the  perturbation  9'  on  element  e,  A9'0  is  A 9'e  at  time  t  =  0.  Note:  the  viscosity  ye 
is  constant  within  each  element  e,  but  is  non-constant  for  different  elements.  For  more 
details  we  refer  to  [30]. 

The  results  for  these  simulations  are  shown  in  Figure  6  for  the  new  DEG  method  a)  and 
for  the  DG  method  with  the  Rusanov  flux  function  b).  Numerical  experiments  clearly 
indicate  that  the  DEG  method  yields  a  more  stable  solution.  Solution  obtained  by  the 
Rusanov  flux  model  exhibits  oscillations  near  the  background  temperature  (9  =  300K), 
which  are  not  present  in  the  simulations  obtained  when  using  the  new  multidimensional 
evolution  Galerkin  operator.  To  get  rid  of  these  oscillations  in  the  Rusanov  model  we 
increased  the  artificial  viscosity  threshold  p0  from  0.1  to  0.7m2/s  (cf.  b)-d)  in  Figure  7). 
Note  that  this  can  influence  large  scale  structures  of  the  solution,  too  (cf.  Test  1). 

In  the  adaptive  viscosity  approach,  the  local  viscosity  on  each  finite  element  scales  ac¬ 
cording  to  (4.18)  and  for  strong  temperature  gradients  it  can  achieve  large  values.  To 
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Figure  6:  Excess  potential  temperature  9'  for  the  Robert  test  [33],  obtained  by  the 
new  semi-implicit  method  and  the  adaptive  artificial  viscosity  according  to  (4.18)  with 
/i0  =  0.1m2/s:  a)  multidimensional  evolution  Galerkin  operator  for  the  numerical  flux 
function  (DEG);  b)  DG  method  with  the  Rusanov  flux  function.  The  real-world  domain 
is  1km x  lkm,  mesh  is  regular,  the  resolution  level  n  =  8,  the  spatial  resolution  fa  22.1m, 
the  shortest  element  edge  ~  44.2m.  The  simulation  times  are  as  indicated.  Note:  for  the 
Rusanov  flux  the  snapshots  are  shown  for  the  last  two  times  only,  where  the  oscillation  of 
the  background  becomes  clearly  noticeable.  Contour  levels  correspond  to  9'  =  —0.05,  0.05, 
0.15,  0.25,  0.35,  and  0.45  K.  Integration  time  step  At  =  0.1s  corresponds  to  CFL  m  1.58, 
advective  CFLU  e  {0,0.0117,0.0123,0.0146}  in  a)  and  CFLU  e  {0.0123,  0.0173}  in  b). 
CFLeg  «  0.0197. 
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Figure  7:  Robert  test  with  the  adaptive  viscosity  according  to  (4.18):  a)  DEG  method 
with  no  =  0.1m2/s;  b)-d)  DG  method  with  the  Rusanov  flux  function  and  po  = 
0.1, 0.4, 0.7m2/s,  respectively.  Mesh  is  regular,  the  resolution  level  n  =  8,  the  spa¬ 
tial  resolution  ps  22.1m,  the  shortest  element  edge  ~  44.2m.  Contour  curves  are  for 
9'  =  —0.05,0.05,0.25,0.4,  CFL  ~  1.57,  advective  CFLU  ~  0.0123  in  a)  and  b), 
CFLa  ps  0.0122  in  c),  and  CFLU  ps  0.0121  in  d).  Simulation  time  t  =  600s. 
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Figure  8:  Robert  test  with  constant  viscosity:  upper  row)  DEG  method  with  fi  —  0  and 
0.1m2/s;  bottom  row  DG  method  with  the  Rusanov  flux  function  and  /i  =  0.2  and  0.4m2/s. 
Mesh  is  regular,  the  resolution  level  n  =  8,  the  spatial  resolution  fa  22.1m,  the  shortest 
element  edge  fa  44.2m.  Contour  curves  are  for  9'  =  —0.05,0.05,0.25,0.4;  CFL  fa  1.57, 
advective  CFLU  fa  0.0132,0.0127,0.0127,0.0122  in  a)-d),  respectively,  CFLeg  ~  0.0197. 
Simulation  time  t  =  600s. 
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Figure  9:  Robert  test  on  the  h— adaptive  mesh  with  adaptive  viscosity  according  to 
(4.18)  with  p0  —  0.1m2/s  for  the  DEG  method  a)  and  the  DG  method  with  the  Ru¬ 
sanov  flux  b).  The  coarse/fine  mesh  resolution  is  n  =  1  —  11,  corresponding  to  the 
spatial  resolution  of  ~  7.8m,  the  shortest  element  edge  ~  15.6m.  Contour  curves  are 
for  9'  =  —0.05,0.05,0.25,0.4;  CFL  &  4.45,  advective  CFLU  &  0.0378  (DEG)  and 
CFLa  0.0379  (Rusanov),  CFLeg  ~  0.055.  Simulation  time  t  =  600s. 


understand  solution’s  behaviour  for  a  certain  given  value  of  the  viscosity  parameter  we 
compare  in  Figure  8  the  DEG  and  the  Rusanov  model  for  fixed  values  of  p.  In  a  fully  in- 
viscid  flow  regime  (when  p  =  0)  and  up  to  /i  =  0.1  the  DEG  model  yields  very  reasonable 
results,  stable  against  the  small  scale  oscillations,  although  some  tiny  island-like  areas  can 
be  observed  around  the  outer  interfaces,  if  compared  with  the  adaptive  artificial  viscosity 
results  in  Figure  7  a).  Hence,  large  scale  structures  can  be  studied  even  for  a  very  low 
artificial  viscosity  and  we  can  analyse  how  small  scale  oscillations  on  the  outer  interface 
are  developing  due  to  the  Kelvin-Helmholtz  instability.  In  the  same  viscosity  regime  the 
results  obtained  by  the  Rusanov  model  are  strongly  oscillating  (not  shown  here)  and  we 
had  to  use  larger  adaptive  viscosity  threshold  parameter  p0  hi  (4.18)  to  obtain  stable 
results,  see  Figures  8  c)  and  d). 

Finally,  in  Figure  9  we  compare  the  DEG  model  and  the  Rusanov  model  using  both 
the  adaptive  viscosity  and  the  adaptive  mesh.  The  adaptive  viscosity  threshold  is  /i0  = 
0.1m2/s  and  the  mesh  resolution  level  changes  between  n  —  1  and  n  =  11.  The  DEG 
method  yields  a  slightly  smoother  solution.  The  outer  interfaces  for  the  Rusanov  model 
exhibit  some  oscillations  on  long-wave  lengths,  but  both  solutions  are  comparable.  The 
tiny  island-like  areas  along  the  outer  interface  on  the  right-hand  side  is  the  trace  of  the 
cold  bubble.  This  experiment  clearly  demonstrates  that  mesh  adaptivity  is  capable  of 
additional  adaptive  numerical  viscosity  and  stabilization  of  a  numerical  solution. 

We  further  analyse  the  convergence  of  the  numerical  scheme  toward  the  exact  solution 
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when  the  grid  resolution  becomes  hirer  by  the  experimental  order  of  convergence  (EOC) 


EOC 
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where  u  is  an  exact  solution,  un  un+\  are  the  solutions  obtained  on  grids  with  the  resolution 
levels  n  and  n  +  1,  respectively.  Since  the  exact  solution  is  unknown  in  our  experiments 
we  take  instead  the  solutions  obtained  on  the  grid  with  higher  resolution  and  calculate  the 
error  between  the  two  solutions  and  EOC  for  the  DEG  and  the  Rusanov  flux  methods. 
They  are  shown  in  Table  1.  The  solution  errors  for  the  DEG  model  are  always  lower 
than  for  the  Rusanov  model,  approximately  by  a  factor  of  1.5  —  2.  The  EOC  of  the  DEG 
model  is  lower  than  in  the  Rusanov  model  on  very  coarse  grids  and  becomes  higher  on 
fine  resolution  grids. 


4.3  Test  3:  Density  current  caused  by  falling  cold  air  bubble 
due  to  Straka  et  al.  [39] 

In  the  density  current  experiment  proposed  by  Straka  et  al.  [39],  the  cold  air  bubble  is 
placed  at  ( xc,yc )  =  (Om,  3000m)  in  the  computational  domain  of  25.6kmx6.4km,  with 
the  maximal  initial  temperature  perturbation  A Tc  =  — 15  K  and  the  distribution  of  the 
initial  temperature  perturbation  according  to 

{0  for  r  >  ry 

(ATc/2)  [1  +  cos  (7rr/rc)]  for  r  <  rc 

The  initial  potential  temperature,  O',  is  calculated  from  AT  and  the  Exner  function,  n, 
using  the  relation  T  =  n6. 

Since  the  density  of  the  cold  air  is  higher,  the  bubble  sinks  gradually  to  the  bottom  of 
the  simulation  domain  (negative  buoyancy)  and  continues  the  viscous  motion  along  the 
bottom  domain  boundary.  Recall  that  we  have  used  the  no-flux  boundary  conditions  in 
our  simulations.  In  the  test  shown  in  Figure  10  we  choose  the  constant  artificial  viscosity 
H  =  75m2/s  from  [39]  and  a  regular  grid  of  the  resolution  level  n  =  8,  since  we  first 
want  to  study  numerical  solutions  without  any  impact  of  the  adaptive  mesh  on  possible 
instabilities.  The  time  integration  was  performed  with  the  time  step  At  =  Is,  which 
corresponds  to  the  following  stability  condition  numbers  CFL  fa  1.365,  CFLU  a  0.141. 
For  this  time  step  r  has  been  chosen  in  such  a  way  that  CFL^q  a  1.7  ■  10~3.  One 
can  see  that  the  new  DEG  method  reproduces  the  flow  in  a  slightly  more  stable  way 
than  the  one-dimensional  Rusanov  flux  model.  At  later  times,  the  oscillations  near  the 
background  temperature  are  stronger  at  this  (and  lower)  resolution  of  the  non- adaptive 
mesh,  even  though  the  viscosity  is  very  high  in  this  test.  Recall  that  the  background 
potential  temperature  has  been  set  to  6  —  300K  in  our  model  and  the  maximal  initial 
potential  perturbation  is  about  9'c  =  — 16. 6  K  in  this  test.  When  we  allow  the  adaptivity 
of  the  mesh  and  choose  a  high  resolution  of  the  grid,  both  models  yield  very  similar 
results,  as  shown  in  Figure  11  for  the  h— adaptive  grid  with  coarse/fine  resolution  levels 
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Figure  10:  Excess  potential  temperature  9'  for  the  density  current  experiment  of  Straka 
et  al.  [39]  obtained  on  a  regular  mesh:  a)  multidimensional  evolution  Galerkin  operator 
for  the  numerical  flux  function  (DEG);  b)  DG  method  with  the  Rusanov  flux  function. 
Artificial  viscosity  is  constant  n  =  75m2/s.  The  real-world  domain  is  25.6kmx6.4km 
(only  part  is  shown),  the  mesh  resolution  level  n  =  8,  the  spatial  resolution  ps  283m, 
the  hnite  element  shortest  edge  ~  566m.  Contour  levels  correspond  to  6'  =  —16  to 
—  1  by  a  step  of  IK.  Note:  the  range  of  the  color  bar  and  colors  in  these  snapshots 
and  in  Figure  11  has  been  restricted  to  — 5K  to  resolve  better  the  fluctuations  near  the 
background  9'  ps  OK  (background  potential  temperature  is  9  =  300K).  Integration  time 
step  At  =  Is,  CFL  ps  1.365,  advective  CFLU  ps  0.141  in  a)  and  CFLU  ps  0.143  in  b), 
CFLeg  ~  1-7  •  10~3.  The  simulation  times  are  as  indicated. 
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Figure  11:  Excess  potential  temperature  9'  for  the  density  current  experiment  of  Straka 
et  al.  [39]  obtained  on  the  adaptive  mesh:  a)  multidimensional  evolution  Galerkin  op¬ 
erator  for  the  numerical  flux  function  (DEG);  b)  DG  method  with  the  Rusanov  flux 
function.  Artificial  viscosity  is  constant  /r  =  75m2/s.  Mesh  coarse/hne  resolution  lev¬ 
els  n  =  1  —  11,  the  finest  spatial  resolution  ss  100m,  the  finite  element  shortest  edge 
ss  200m.  Contour  levels  correspond  to  6'  =  —16  to  —1  by  a  step  of  IK.  Integration 
time  step  At  =  0.5s,  CFL  ~  1.91,  advective  CFLU  e  {0,0.194,0.182,0.179}  in  a)  and 
CFLU  e  {0,0.194,0.182,0.178}  in  b).  CFLeg  ~  2.4  •  10~3.  The  simulation  times  are  as 
indicated. 
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a)  BDF2+EG: 


time  =  100s 


n  =  gridlevel 

|  \un—  un+2 1|  /volume 

||un+2  -  un+4  II  /volume 

EOC(n,  n  +  2,  n  +  4) 

3 

0.056854 

0.028865 

0.9779 

4 

0.046331 

0.019119 

1.2769 

5 

0.028865 

0.006807 

2.0842 

6 

0.019119 

0.002648 

2.8519 

7 

0.006807 

0.000866 

2.9740 

time  =  150s 


n  =  gridlevel 

||  un-  un+2 1|  /volume 

||un+2  -  un+4  II  /volume 

EOC{n,  n  +  2,  n  +  4) 

3 

0.092453 

0.047843 

0.9504 

4 

0.075194 

0.030815 

1.2870 

5 

0.047843 

0.011797 

2.0199 

6 

0.030815 

0.004206 

2.8730 

7 

0.011797 

0.001539 

2.9383 

b)  BDF2+Rusanov: 


time  =  100s 


n  =  gridlevel 

| \un-  un+2 1| /volume 

||un+2  -  un+4 1|  /volume 

EOC(n,  n  +  2,  n  +  4) 

3 

0.068548 

0.033892 

1.0162 

4 

0.059404 

0.023547 

1.3350 

5 

0.033892 

0.010428 

1.7004 

6 

0.023547 

0.004046 

2.5522 

7 

0.010428 

0.001697 

2.6195 

time  =  150s 


n  =  gridlevel 

||  un~  un+2 1|  /volume 

\\un+2  ~  Un+4 1| /volume 

EOC(n ,  n  +  2,  n  +  4) 

3 

0.126310 

0.065635 

0.9445 

4 

0.107820 

0.040578 

1.4099 

5 

0.065635 

0.020024 

1.7128 

6 

0.040578 

0.007314 

2.4719 

7 

0.020024 

0.003247 

2.6246 

Table  1:  The  solution  error  (Z/2-norm)  and  the  experimental  order  of  convergence  (EOC) 
in  the  large  time  step  simulations  with  a)  BDF2+EG  and  b)  BDF2+Rusanov  of  the 
Robert  test  [33]  with  constant  viscosity  /i  =  0.1m2/s. 


n  =  1  —  11.  Here  both  the  DEG  and  Rusanov  flux  models  perform  very  well,  the  differences 
between  the  solutions  are  hardly  distinguishable  anymore.  In  this  simulation  the  time  step 
At  =  0.5s  was  reduced  by  a  factor  of  two  to  adjust  it  to  the  finer  grid,  if  compared  to  the 
experiment  with  the  regular  mesh  n  =  8. 
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5  Conclusions 


In  the  present  paper  we  have  derived  a  new  adaptive  discontinuous  evolution  Galerkin 
method.  The  novelty  of  our  approach  relies  in  the  combination  of  the  discontinuous 
Galerkin  method  with  a  genuinely  multidimensional  evolution  operators  based  on  the 
theory  of  bicharacteristics  for  underlying  hyperbolic  balance  laws.  In  this  paper  the  DEG 
method  is  applied  for  test  cases  modeling  dry  atmospheric  convection.  In  order  to  take 
into  account  multiscale  phenomena  that  typically  arise  in  atmospheric  flows  we  split  fluxes 
into  a  linear  part  governing  acoustic  waves  and  the  resulting  nonlinear  part.  The  linear 
operator  has  to  be  chosen  in  such  a  way  that  the  fastest  waves  of  the  system  are  retained, 
although  in  their  linearized  form.  Time  integration  is  realized  by  the  IMEX  type  approx¬ 
imation  using  the  semi-implicit  BDF2  method.  In  order  to  efficiently  resolve  small  scale 
flow  structures  adaptive  mesh  refinement  is  used.  This  is  realized  via  the  AMATOS  func¬ 
tion  library.  Numerical  experiments  presented  in  Section  4  demonstrate  high  accuracy, 
stability  and  robustness  of  the  new  method  and  illustrate  that  complex  multidimensional 
flow  structures  are  approximated  in  a  better  way  than  by  the  discontinuous  Galerkin 
method  with  a  standard  one-dimensional  numerical  flux,  e.g.,  the  Rusanov  flux  function. 
The  Rusanov  flux  model  is  fast  and  can  yield  solutions  of  similar  quality  as  the  DEG 
model  when  adaptive  mesh  refinement  and  adaptive  artificial  viscosity  are  used.  On  the 
other  hand,  the  DEG  method  is  more  stable  due  to  the  truly  multidimensional  nature  of 
the  EG  operator.  In  practice,  this  implies  that  less  or  no  artificial  viscosity  is  required.  For 
low  viscosity  regimes  and  for  coarse  grids,  where  the  Rusanov  flux  model  can  be  unstable, 
the  DEG  model  performs  better.  Further  increase  of  the  efficiency  of  the  DEG  method 
can  be  achieved  by  porting  the  calculation  of  the  EG  operator  to  graphics  processing  units 
(GPU),  as  has  been  done  recently  for  the  explicit  DEG  method  [6].  In  the  future  it  will 
be  interesting  to  generalize  the  DEG  method  for  three-dimensional  flows  and  apply  it  to 
more  complex  atmospheric  problems,  such  as  simulation  of  a  cloud  environment. 
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